package scu.maqiang.heat;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES3T41;
import scu.maqiang.mesh.Mesh3T4;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.SRMatrix;

import java.util.Arrays;

public class Heat3DT4 {
    public static void main(String[] args) {
        Mesh3T4 mesh = new Mesh3T4();
        mesh.cube(20, 20, 20);

        FES3T41 fs = new FES3T41(mesh);
        int dof = fs.GetNdof();
        SRMatrix A = new SRMatrix(fs.GetNdof());
        SRMatrix M = new SRMatrix(fs.GetNdof());
        fs.assembleStiff(new double[]{1.0}, BVPType.COMMON, A);
        fs.assembleMass(new double[]{1.0}, BVPType.CONSIST_MASS, M);

        double[] RHS = new double[dof];
        double[] SOURCE = new double[dof];
        fs.assembleSource(new double[]{1.0}, BVPType.COMMON, SOURCE);

        double theta = 0.5;
        double T = 1.0;
        double dt = T / 100.0;

        SRMatrix Lmat = M.add(theta * dt, A);
        SRMatrix Rmat = M.add((theta - 1) * dt, A);

        fs.applyBC_MBN(Lmat, RHS, Direct.X, 0.0, 1, 2, 3, 4, 5, 6);
        NewIterSSolver solver = new NewIterSSolver(Lmat);

        double time = 0.0;
        double[] x = new double[dof];
        mesh.toTecplot("Heat.dat", time, x);
        while(time < T) {
            time += dt;
            double[] temp = Rmat.mul(x);
            Arrays.setAll(RHS, i-> SOURCE[i] * dt + temp[i]);
            fs.applyBC_MBN(RHS, Direct.X, 0.0, 1, 2, 3, 4, 5, 6);
            solver.PCGSSOR(RHS, x, 1.5, 1);
            mesh.toTecplot("Heat.dat", time, x);
        }
    }

}
